clear
set more off


*Get data
use ${tempdir}/panel_trends_30yr_income.dta, clear
keep if race==1


*Variables for the relevant marriagemarket and childhood years
gen maryear = year2 - 10
gen chyear = year2 - 30

*Here, assign the state the match the panel controls as statefip, state of residence
*Change to bpl to match on birth place instead
gen fips = statefip

* Regressions of ever married status
egen firstcode = group(sex first year2)

gen int bpl2 = bpl
replace bpl2 = 9999 if bpl>100

gen byte foreign_born = bpl2>100 if bpl2~=.

sort year2 foreign_born

gen logoccscore_gap = logoccscore_father_spouse - logoccscore_father
gen abs_agegap = abs(agegap)

gen urb = urban==2 if urban!=.
label var urb "Urban dummy"


foreach y in $yearlist {
		gen y`y' = year2==`y'
}


local dep_ever_married = 1
local dep_logocc_fil = 2


gen logocc_fil = logoccscore_father_spouse 
gen logoccscore_husband_temp = log(husboccscore)
gen husb_parent_qtl = qtl_score_wgt_full_spouse

local qs4 $qs4
local qs6 $qs6
local qs8 $qs8
local qs10 $qs10 

rename qtl_score_wgt_full qtl4_score_wgt_full

cap postclose results
postfile results tag nq depvar spec qtile year b se  using "$tempdir/reg_robust_differentquantiles.dta", replace


foreach var of varlist ever_married logocc_fil  {

	foreach qm in 4 6 8 10 {
	
	cap drop q1* q2* q3* q4*
	cap drop q5* q6*
	cap drop q7* q8*
	cap drop q9* 
	cap drop q10*
	
	di "beginning of `qm'"
	
	forvalues i = 1/`qm' {
	gen q`i' = qtl`qm'_score_wgt_full==`i'
	replace q`i'=. if qtl`qm'_score_wgt_full==.
}

forvalues i=1/`qm' {
	foreach y in $yearlist {
		gen q`i'_`y' = q`i'*y`y'
	}
}

* Ever married regressions -- no bpl controls
reg `var' `qs`qm''  i.age foreign_born if sex==2, cluster(firstcode)
gen s1_`var'1_q`qm' = e(sample)

forvalues i=2/`qm' {
	foreach y in $yearlist {
		local b`i'_`y' = _b[q`i'_`y']
		local se`i'_`y' = _se[q`i'_`y']
	}
}

forvalues i=2/`qm' {
	foreach y in $yearlist {
		post results (1) (`qm') (`dep_`var'') (1) (`i') (`y') (`b`i'_`y'') (`se`i'_`y'') 
	}
}


* Ever married regressions -- full bpl controls
reg `var' `qs`qm'' i.age i.bpl2##i.year2 if sex==2, cluster(firstcode)
gen s2_`var'1_q`qm' = e(sample)



forvalues i=2/`qm' {

	foreach y in $yearlist{
		local b`i'_`y' = _b[q`i'_`y']
		local se`i'_`y' = _se[q`i'_`y']
	}
}

forvalues i=2/`qm' {
	foreach y in $yearlist {
		post results (1) (`qm') (`dep_`var'') (2) (`i') (`y') (`b`i'_`y'') (`se`i'_`y'') 
}

}
}
}

postclose results

save "${tempdir}reg_data_quantiles.dta", replace 

clear 
save "${tempdir}reg_nobs_quantiles.dta", replace empty

local dep_ever_married = 1
local dep_logocc_fil = 2


foreach x in 1 {
foreach var in ever_married logocc_fil  {
foreach spec in 1 2 {
foreach qm in 4 6 8 10 {

use "${tempdir}reg_data_quantiles.dta", clear

keep if s`spec'_`var'`x'_q`qm' == 1

sort year2 first
by year2 first: gen nnames_female = _n == 1

if "`var'" == "logocc_fil" {
sort year2 spousefirst
by year2 spousefirst: gen nnames_male = _n == 1
	}
	
collapse (sum) nobs=s`spec'_`var'`x'_q`qm' nnames_*, by(year2)	

gen tag = `x'
gen depvar = `dep_`var''
gen spec = `spec'
gen nq = `qm'

rename year2 year 

append using "${tempdir}reg_nobs_quantiles.dta"
save "${tempdir}reg_nobs_quantiles.dta", replace 

}
}
}
}

use "$tempdir/reg_robust_differentquantiles.dta", clear
merge m:1 tag depvar spec nq year using "${tempdir}reg_nobs_quantiles.dta"

export excel using "$outputdir/Regression Robustness.xlsx", sheet("differentquantiles", replace) firstrow(var)

erase "${tempdir}reg_data_quantiles.dta"
erase "${tempdir}reg_nobs_quantiles.dta"

